In search of the Goldilocks zone for hybrid speciation II: hard times for hybrid speciation?
This is the analysis file corresponding to the manuscript “In search of the Goldilocks zone for hybrid speciation II: hard times for hybrid speciation”.
Reading the simulations output and generating the main dataset files
Default values to debug the add_dataframe() function.
if (FALSE){
k="output_sym_dom_rfunction_pi_0d5_rdot20.txt"
root="output_sym_adjabab_pi_0d5_nos"
}Loading functions and packages
source("functions.R")
library(ggplot2)## Warning: le package 'ggplot2' a été compilé avec la version R 4.2.2print(paste(c("ggplot version ", as.character(packageVersion("ggplot2"))),collapse=""))## [1] "ggplot version 3.4.1"library(fields)## Warning: le package 'fields' a été compilé avec la version R 4.2.2## Le chargement a nécessité le package : spam## Warning: le package 'spam' a été compilé avec la version R 4.2.2## Spam version 2.9-1 (2022-08-07) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attachement du package : 'spam'## Les objets suivants sont masqués depuis 'package:base':
##
## backsolve, forwardsolve## Le chargement a nécessité le package : viridis## Warning: le package 'viridis' a été compilé avec la version R 4.2.2## Le chargement a nécessité le package : viridisLite## Warning: le package 'viridisLite' a été compilé avec la version R 4.2.2##
## Try help(fields) to get started.print(paste(c("fields version ", as.character(packageVersion("fields"))),collapse=""))## [1] "fields version 14.1"library(viridis)
print(paste(c("viridis version ", as.character(packageVersion("viridis"))),collapse=""))## [1] "viridis version 0.6.2"library(RColorBrewer)
print(paste(c("RColorBrewer version ", as.character(packageVersion("RColorBrewer"))),collapse=""))## [1] "RColorBrewer version 1.1.3"library(cowplot)## Warning: le package 'cowplot' a été compilé avec la version R 4.2.2print(paste(c("cowplot version ", as.character(packageVersion("cowplot"))),collapse=""))## [1] "cowplot version 1.1.1"library(plyr)## Warning: le package 'plyr' a été compilé avec la version R 4.2.2print(paste(c("plyr version ", as.character(packageVersion("plyr"))),collapse=""))## [1] "plyr version 1.8.8"Load default dataset (1,000 replicates per combination of parameter)- Per default load the resulting dataset directly.
if (FALSE){
whole_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
colnames(whole_dataset)=names_col
known_folder=c()
} else{
whole_dataset=read.csv("whole_dataset.csv")
#whole_dataset$arch=as.factor(list_Arch[whole_dataset$arch])
known_folder=read.table("know_folder.txt")
}Retrieve all possible folders starting by “output_” in the working directory.
list_folder=list.dirs(path = ".", full.names = TRUE, recursive = TRUE)
list_folder=gsub("./","",list_folder)
list_folder=list_folder[grep("^output", list_folder)]Check the current main dataframe and update it if necessary
update=FALSE
for (i in list_folder){
if (!(i %in% unlist(known_folder))){
update=TRUE
#check_non_finish(i)
temp_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
temp_dataset2=add_dataframe(i,temp_dataset,truncature = 1000)
temp_dataset2=temp_dataset2[-c(1,2),]
colnames(temp_dataset2)=names_col
if(!check_arch(temp_dataset2)){
print(i)
break
}
whole_dataset=rbind(whole_dataset,temp_dataset2)
known_folder=c(known_folder,i)
}
if(length(known_folder)%%40==0){print(length(known_folder)/length(list_folder))}
}
whole_dataset=whole_dataset[!is.na(whole_dataset$n_rep),]
if (update){
write.table(known_folder,"know_folder.txt")
write.csv(whole_dataset,"whole_dataset.csv",row.names = FALSE)
}Reload the data for better consistency
whole_dataset=read.csv("whole_dataset.csv")
whole_dataset$arch=as.factor(whole_dataset$arch)
whole_dataset=whole_dataset[-c(1,2),]Load detailed dataset (10,000 replicates per combination of parameter) Per default load the resulting dataset directly.
if (FALSE){
whole_dataset_extra=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
colnames(whole_dataset_extra)=names_col
known_folder_extra=c()
dismiss_folder_extra=c()
} else{
whole_dataset_extra=read.csv("whole_dataset_extra.csv")
#whole_dataset_extra$arch=as.factor(list_Arch[whole_dataset_extra$arch])
known_folder_extra=read.table("know_folder_extra.txt")
}Check the current main dataframe and update it if necessary
update=FALSE
temp_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
for (i in list_folder){
if (!(i %in% unlist(known_folder_extra))){
update=TRUE
temp_dataset2=add_dataframe(i,temp_dataset,truncature = 10000)
temp_dataset2=temp_dataset2[-c(1,2),]
colnames(temp_dataset2)=names_col
if(!check_arch(temp_dataset2)){
print(i)
break
}
whole_dataset_extra=rbind(whole_dataset_extra,temp_dataset2)
known_folder_extra=c(known_folder_extra,i)
}
if(length(known_folder_extra)%%40==0){print(length(known_folder_extra)/length(list_folder))}
}
whole_dataset_extra=whole_dataset_extra[!is.na(whole_dataset_extra$n_rep),]
if (update){
write.table(known_folder_extra,"know_folder_extra.txt")
write.csv(whole_dataset_extra,"whole_dataset_extra.csv",row.names = FALSE)
}Reload the data for better consistency
whole_dataset_extra=read.csv("whole_dataset_extra.csv")
whole_dataset_extra$arch=as.factor(whole_dataset_extra$arch)
whole_dataset_extra=whole_dataset_extra[-c(1,2),]Analysis of the unlinked loci case
General outcomes
Relationship between survival, epistasis and initial population size. (Since all architectures are equivalent, we always used the Adjacent ABAB architecture for the case of unlinked loci.)
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB","sin_DMI")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0&whole_dataset$e12<0,]
temp_data$err_surv=sqrt((temp_data$p_alive/temp_data$n_rep)*(1-(temp_data$p_alive/temp_data$n_rep))/temp_data$n_rep)#pdf("Figure/P_alive_vs_ep_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
temp_name2=c(`0`="One DMI",`1`="Two DMIs")
names(temp_name2)=c("0","1")
extra<-data.frame(label=c("A","B","C","D"),h12=c(0,1,0,1),arch=c("sin_DMI","sin_DMI","adj_ABAB","adj_ABAB"))
# adding extra column
temp_data$Cond=sapply(c(1:length(temp_data[,1])), function(x) paste(temp_data$h12[x], temp_data$arch[x], sep="_"))
temp_data$extra=plyr::mapvalues(temp_data$Cond, c("0_sin_DMI", "1_sin_DMI" , "0_adj_ABAB","1_adj_ABAB"), c("A","B","C","D"))
#str(temp_data)
ggplot(temp_data,aes(x=-e12,y=p_alive/n_rep,col=as.factor(N)))+geom_point()+scale_x_log10()+ xlab("Strength of epistasis")+ylab("Prob.(Survival)")+geom_line(size=0.5)+labs(col = "Initial pop. size")+geom_text(x=-4.3,y=0.9,aes(label=extra), colour="black")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name),)+theme(legend.position="bottom")+ylim(0,1)+theme(text = element_text(size=14))+theme(panel.spacing.x = unit(1, "lines"))+
coord_cartesian(clip = "off")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
#dev.off()temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("sin_DMI","adj_ABAB")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0&whole_dataset$e12<0&whole_dataset$h12==1,]
temp_data$f1_fitness=(1+temp_data$e12)*(1+temp_data$e34)
temp_data$arch=as.character(temp_data$arch)
temp_data$arch[temp_data$arch=="sin_DMI"]="Single DMI"
temp_data$arch[temp_data$arch=="adj_ABAB"]="Two DMIs"
# adding extra column
temp_data$Cond=sapply(c(1:length(temp_data[,1])), function(x) paste(temp_data$h12[x], temp_data$arch[x], sep="_"))
temp_data$extra=plyr::mapvalues(temp_data$Cond, c("0_sin_DMI", "1_sin_DMI" , "0_adj_ABAB","1_adj_ABAB"), c("A","B","C","D"))## The following `from` values were not present in `x`: 0_sin_DMI, 1_sin_DMI, 0_adj_ABAB, 1_adj_ABAB#str(temp_data)
#pdf("Figure/one_vs_two_DMI_dame_hybrid_load.pdf",width=8, height=6)
ggplot(temp_data,aes(x=1-f1_fitness,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+ xlab("Hybrid load of the F1 individuals")+ylab("Prob.(Survival)")+geom_line(size=0.5)+labs(col = "Number of DMIs")+facet_wrap(temp_data$N,nrow=2)+theme(legend.position="bottom")+ylim(0,1)+xlim(0.01,1)+theme(panel.spacing.x = unit(1, "lines"))## Warning: Removed 64 rows containing missing values (`geom_point()`).## Warning: Removed 8 rows containing missing values (`geom_line()`).
#dev.off()Effect of epistasis on the resolution time of the genetic conflict of a single DMI and two DMIs.
#pdf("Figure/Res_time_vs_ep_unlinked.pdf",width=8, height=6)
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("sin_DMI")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0,]
ggplot(temp_data,aes(x=-e12,y=t_alive_dmi1,col=as.factor(N)))+geom_point()+scale_x_log10()+ xlab("Strength of epistasis")+ylab("Resolution time")+labs(col = "Initial pop. size")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+geom_line()+scale_y_log10()+theme(text = element_text(size=14))+theme(panel.spacing.x = unit(1, "lines"))## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$K==1000000&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000)&whole_dataset$s1<0,]
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
temp_name2=c(`0`="One DMI",`1`="Two DMIs")
names(temp_name2)=c("0","1")
ggplot(temp_data,aes(x=-e12,y=t_alive_both,col=as.factor(N)))+geom_point()+scale_x_log10()+ xlab("Strength of epistasis")+ylab("Resolution time")+labs(col = "Initial pop. size")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+geom_line()+scale_y_log10()+theme(text = element_text(size=14))+theme(panel.spacing.x = unit(1, "lines"))## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
#geom_smooth(se=F)+ xlim(-1,0)+
#dev.off()Relationship between survival, epistasis and initial population size, as above, but with population size on the X axis.
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$K==1000000&whole_dataset$s1==-0.001&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB","sin_DMI")&whole_dataset$e12 %in% c(0,-0.001,-0.01,-0.1,-0.5,-0.99),]#pdf("Figure/P_alive_vs_Nini_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
temp_name2=c(`0`="One DMI",`1`="Two DMIs")
names(temp_name2)=c("0","1")
extra<-data.frame(label=c("A","B","C","D"),h12=c(0,1,0,1),arch=c("sin_DMI","sin_DMI","adj_ABAB","adj_ABAB"))
# adding extra column
temp_data$Cond=sapply(c(1:length(temp_data[,1])), function(x) paste(temp_data$h12[x], temp_data$arch[x], sep="_"))
temp_data$extra=plyr::mapvalues(temp_data$Cond, c("0_sin_DMI", "1_sin_DMI" , "0_adj_ABAB","1_adj_ABAB"), c("A","B","C","D"))
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(-e12)))+geom_point(size=1)+geom_line()+scale_x_log10()+ xlab("Initial pop. size")+ylab("Prob.(Survival)")+labs(col = "Strength of epistasis")+facet_grid(as.numeric(temp_data$arch!="sin_DMI")~h12,labeller =labeller(.rows=temp_name2,.cols=temp_name))+theme(legend.position="bottom")+ylim(0,1)+geom_text(x=2,y=.97,aes(label=extra), colour="black")+theme(text = element_text(size=14))+theme(panel.spacing.x = unit(1, "lines"))
#dev.off()Relationship between minimum population size and epistasis
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000),]#pdf("Figure/Min_pop_size_vs_ep_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
ggplot(temp_data,aes(x=-e12,y=N_min_alive,col=as.factor(N)))+geom_point()+geom_line()+ xlab("Strength of epistasis")+ylab("Minimum pop. size")+labs(col = "Initial pop. size")+xlim(0,1)+facet_wrap(~h12,labeller =as_labeller( temp_name))+theme(legend.position="bottom")+coord_cartesian(ylim = c(0, 500))+scale_x_log10()## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.## Warning: Transformation introduced infinite values in continuous x-axis
## Transformation introduced infinite values in continuous x-axis
#dev.off()Relationship between resolution time and epistasis
#pdf("Figure/Time_resolution_vs_ep_unlinked.pdf",width=8, height=6)
temp_name=c(`0`="Recessive",`1`="Codominant")
names(temp_name)=c("0","1")
ggplot(temp_data,aes(x=-e12,y=t_both,col=as.factor(N)))+geom_point() +geom_line()+ xlab("Strength of epistasis")+ylab("Time resolution 2 DMIs")+labs(col = "Initial pop. size")+xlim(0,1)+facet_wrap(~h12,labeller =as_labeller( temp_name))+theme(legend.position="bottom")+coord_cartesian(ylim = c(0, 100))
#dev.off()temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$h12==1 &whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000),] # codominant case
temp_data2=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$h12==0 &whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000),]#recessive caseggplot(temp_data,aes(x=N_min_alive,y=t_both,col=p_alive/n_rep))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 codominant DMIs")+labs(col = "P(alive)")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
ggplot(temp_data,aes(x=N_min_alive,y=t_both,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 codominant DMIs")+labs(col = "Initial pop. size")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
ggplot(temp_data2,aes(x=N_min_alive,y=t_both,col=p_alive/n_rep))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 recessive DMIs")+labs(col = "P(alive)")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
ggplot(temp_data2,aes(x=N_min_alive,y=t_both,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution 2 recessive DMIs")+labs(col = "Initial pop. size")+theme(legend.position="bottom")+scale_x_log10()+scale_y_log10()
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$r12==0.5&whole_dataset$h12==1 &whole_dataset$arch %in% c("adj_ABAB")&whole_dataset$N %in% c(100,500,1000,2500,5000,10000)&whole_dataset$e12< -0.09,] # codominant case#pdf("Figure/N_min_vs_t_res_unlinked.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N_min_alive,y=t_both,col=as.factor(e12)))+geom_point()+ xlab("Minimum population size")+ylab("Time to resolution of 2 DMIs")+labs(col = "Initial pop. size")+theme(legend.position="bottom",text = element_text(size=14))+scale_x_log10()+scale_y_log10()+facet_wrap(~N)+labs(tag = "A")
#dev.off()Deterministic dynamics
To better understand the outcome observed above, and why intermediate epistasis values presented the highest risk of extinction, we focused on the deterministic scenario. ### For an initial pop size of N=5000
data0=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep0.txt")
data1=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.1.txt")
data2=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.2.txt")
data3=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.3.txt")
data4=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.4.txt")
data5=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.5.txt")
data6=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.6.txt")
data7=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.7.txt")
data8=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.8.txt")
data9=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.9.txt")
data99=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N5000/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N5000_ep-0.99.txt")Population size evoltuion for different epistasis strengths
#pdf("Figure/Determ_pop_size_N5000_unlinked.pdf",height=4.5,width=4)
plot(data5$V1,type='l',ylim=c(3,500),xlim=c(0,750),col="purple",xlab="Generations",ylab="Pop. size",lwd=2,log="y",cex=2)
points(seq(1000),data6$V1,type='l',col="blue",lwd=2)
points(seq(1000),data7$V1,type='l',col="cyan",lwd=2)
points(seq(1000),data8$V1,type='l',col="orange",lwd=2)
points(seq(1000),data9$V1,type='l',col="red",lwd=2)
legend("bottomright",legend=c(expression(paste(epsilon,"=-0.5")),expression(paste(epsilon,"=-0.6")),expression(paste(epsilon,"=-0.7")),expression(paste(epsilon,"=-0.8")),expression(paste(epsilon,"=-0.9"))),col=c("purple","blue","cyan","orange","red"),pch=16)
text(x=-180, y=740, substitute(paste(bold('C'))), xpd=NA)
#dev.off()We want to calculate here the initial mean growth rate of the population to compare it to en evolutionary rescue scenario (which selection coefficient in a single locus model will generate the same initial and final growth rate).
#pdf("Figure/rbar-vs-ep.pdf",width=8,height=6)
plot(-c(0,0.1,0.2,0.3,0.4,.5,.6,.7,.8,.9,.99),c(data0$V1[2]/data0$V1[1],data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1]),xlab=expression(paste("Epistasis, ",epsilon)),ylab="Average r",pch=16,col="blue",type="b",cex.axis=1.2,cex.lab=1.2,ylim=c(0.4,1.01))
legend("topleft",legend = c("F1","F2"),col=c("blue","cyan"),pch=16,cex=1.2)
points(-c(0,0.1,0.2,0.3,0.4,.5,.6,.7,.8,.9,.99),c(data0$V1[2]/data0$V1[1],data1$V1[3]/data1$V1[2],data2$V1[3]/data2$V1[2],data3$V1[3]/data3$V1[2],data4$V1[3]/data4$V1[2],data5$V1[3]/data5$V1[2],data6$V1[3]/data6$V1[2],data7$V1[3]/data7$V1[2],data8$V1[3]/data8$V1[2],data9$V1[3]/data9$V1[2],data99$V1[3]/data99$V1[2]),col="cyan",pch=16,type="b")
#dev.off()c(data0$V1[2]/data0$V1[1],data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1])## [1] 1.0059661 0.9103993 0.8248922 0.7494448 0.6840570 0.6287288 0.5834604
## [8] 0.5482516 0.5231024 0.5080129 0.5030334sqrt(c(data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1])/1.005966)-1## [1] -0.04868509 -0.09446145 -0.13686613 -0.17537882 -0.20943053 -0.23842262
## [7] -0.26175878 -0.27888967 -0.28936640 -0.29285779#pdf("Figure/equivalent_s.pdf",width=8,height=6)
plot(-c(0,0.1,0.2,0.3,0.4,.5,.6,.7,.8,.9,.99),sqrt(c(1,data1$V1[2]/data1$V1[1],data2$V1[2]/data2$V1[1],data3$V1[2]/data3$V1[1],data4$V1[2]/data4$V1[1],data5$V1[2]/data5$V1[1],data6$V1[2]/data6$V1[1],data7$V1[2]/data7$V1[1],data8$V1[2]/data8$V1[1],data9$V1[2]/data9$V1[1],data99$V1[2]/data99$V1[1])/1.005966)-1,pch=16,xlab="Epistasis",ylab="Equivalent s")
#dev.off()For an intial population size of N=500
Same as above but with a reduction by a factor 10 of the initial population size.
data3=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.3.txt")
data4=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.4.txt")
data5=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.5.txt")
data6=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.6.txt")
data7=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.7.txt")
data8=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.8.txt")
data9=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.9.txt")
data99=read.table("deterministic_output/output_sym_cod_r0d5_pi_0d5_adjabab_N500/freq_output_sym_cod_r0d5_pi_0d5_adjabab_N500_ep-0.99.txt")#pdf("Figure/Determ_pop_size_N500_unlinked.pdf",height=4.5,width=4)
plot(data5$V1,type='l',ylim=c(3,500),xlim=c(0,750),col="purple",xlab="Generations",ylab="Pop. size",lwd=2,log="y",cex=2)
points(seq(1000),data6$V1,type='l',col="blue",lwd=2)
points(seq(1000),data7$V1,type='l',col="cyan",lwd=2)
points(seq(1000),data4$V1,type='l',col="darkorchid4",lwd=2)
points(seq(1000),data3$V1,type='l',col="black",lwd=2)
legend("bottomright",legend=c(expression(paste(epsilon,"=-0.3")),expression(paste(epsilon,"=-0.4")),expression(paste(epsilon,"=-0.5")),expression(paste(epsilon,"=-0.6")),expression(paste(epsilon,"=-0.7"))),col=c("black","darkorchid4","purple","blue","cyan","orange"),pch=16)
text(x=-180, y=740, substitute(paste(bold('B'))), xpd=NA)
#dev.off()#Analysis of the linked loci: Codominant scenario
All 6 possible architectures
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$h12==1&whole_dataset$arch!="sin_DMI",]Legned for the 6 architectures
#pdf("Figure/legend_arch.pdf",width=7,height=1)
plot_grid(get_legend(ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point(size=4)+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=20),legend.position = "bottom")+ylim(0,1)))
#dev.off()Relationship between survival probability and minimum population size for all 6 architectures
#pdf("Figure/Min_pop_size_vs_alive_cod.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob.(Survival)")+labs(col = "Recombination rate")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=t_both,y=p_alive/n_rep,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the 2 DMIs")+ylab("Prob.(Survival)")+labs(col = "Recombination rate")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=log10(t_both)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob.(Survival)")+labs(col = "(log10) Time of resolution of the 2 DMIs")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=N,y=N_min_alive,col=p_alive/n_rep))+geom_point()+geom_jitter(width=0.2)+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Starting population size")+ylab("Minimum population size")+labs(col = "Prob.(Survival)")+scale_x_log10()+scale_y_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom")
ggplot(temp_data,aes(x=t_both,y=N_min_alive,col=p_alive/n_rep))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the 2 DMIs")+ylab("Minimum population size")+labs(col = "Prob.(Survival)")+scale_x_log10()+scale_y_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
ggplot(temp_data,aes(x=t_both,y=p_alive/n_rep,col=log10(N_min_alive)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the 2 DMIs")+ylab("Prob.(Survival)")+labs(col = "Minimum population size")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))
#dev.off()Same as above but with guidelines
#pdf("Figure/Min_pop_size_vs_alive_cod_guide.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=log10(t_both)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob.(Survival)")+labs(col = "(log10) Time of resolution of the 2 DMIs")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))+ geom_segment(aes(x = 10, y = 0, xend = 1000, yend = 1),size=1,color="black")
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=(r12)))+geom_point(size=0.5)+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("Prob.(Survival)")+labs(col = "Recombination rate")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label_long))+ theme(legend.position="bottom",legend.key.size = unit(0.33, "in"))+ geom_segment(aes(x = 10, y = 0, xend = 600, yend = 1),size=1,color="black")+ theme(axis.text=element_text(size=10),axis.title =element_text(size=12) )+theme(strip.text.x = element_text(size = 12))+guides(color=FALSE)## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
#dev.off()temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$h12==1&whole_dataset$arch!="sin_DMI"&whole_dataset$N %in% c(100,500,1000,2500,5000,10000,50000,100000),]#pdf("Figure/p_alive_vs_Nmin.pdf",width=8, height=6)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Initial Population size")+scale_x_log10()+ theme(legend.position="bottom")+scale_colour_brewer(palette="Set1")+ylim(0,1)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Initial Population size")+scale_x_log10()+ theme(legend.position="bottom")+scale_colour_brewer(palette="Set1")+ylim(0,1)
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(N)))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Initial Population size")+scale_x_log10()+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")+scale_colour_brewer(palette="Set1")+ylim(0,1)
#dev.off()Epistasis is fixed to e=-0.2
Relationship between hybrid speciation and recombination for all 6 architectures and different population sizes
#pdf("all_6arch_hard_sel.pdf",width=8,height=6)
#adj ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Adjacent ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#adj ABBA
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Adjacent ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Crosed ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Crossed AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#nes ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Nested ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#nes AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.",pch=16,xlim=c(0.001,0.5), main="Nested AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$n_rep,col="orchid1",pch=16)
#dev.off()Relationship between hybrid speciation conditioned on survival and recombination for all 6 architectures and different population sizes
#pdf("all_6arch_hard_sel_cond.pdf",width=8,height=6)
#adj ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Adjacent ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#adj ABBA
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Adjacent ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="adj_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Crosed ABBA DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_ABBA",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#cro AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Crossed AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="cro_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#nes ABAB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Nested ABAB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_ABAB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#nes AABB
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
plot(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid4",log="x",ylim=c(0,1),xlab="Recombination",ylab="Hybrid speciation pro.|survival",pch=16,xlim=c(0.001,0.5), main="Nested AABB DMIs",cex.axis=1.2,cex.lab=1.2)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==5000&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="darkorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==500&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="mediumorchid2",pch=16)
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$mean_off==1.01&whole_dataset$N==50&whole_dataset$arch=="nes_AABB",]
temp_data=temp_data[!duplicated(temp_data$r12),]
points(temp_data$r12,(temp_data$p_alive_hybA1+temp_data$p_alive_hybB1)/temp_data$p_alive,col="orchid1",pch=16)
#dev.off()Initial population size of N=5000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]Relationship between survival probability, hybrid speciation and recombination for all architectures
#pdf("Figure/cod_case_ep0d2_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylim(0.4,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.25) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid sp.|survival)")+xlab("Recombination")+geom_line()#+ guides(col=FALSE)## Warning: Removed 3 rows containing missing values (`geom_point()`).
#dev.off()relationship between minimum population size, hybrid speciation probbaility and recombination
#pdf("Figure/Nminvsrec_and_hyybvsrec_cod0d2.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
#dev.off()Further exploration of this scenario, including time of resolution of the DMIs
ggplot(temp_data,aes(y=t_alive_both,x=N_min_alive,col=p_alive/n_rep))+geom_point()+ xlab("Minimum population size")+ylab("Time resolution")+labs(col = "P(alive)")+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")+scale_colour_gradient(low = "orange", high = "purple")
ggplot(temp_data,aes(x=N_alive_both,y=t_alive_both,col=p_alive/n_rep))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Time of resolution of the two DMIs")+labs(col = "Prob(surviving)")
ggplot(temp_data,aes(x=N_alive_both,y=t_alive_both,col=((p_hybA1+p_hybB1)/p_alive)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Time of resolution of the two DMIs")+labs(col = "Prob(hyb_spec|alive)")
ggplot(temp_data,aes(x=N_alive_both,y=t_alive_both,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Time of resolution of the two DMIs")+labs(col = "Recombination rate")
ggplot(temp_data,aes(x=N_alive_both,y=p_alive/n_rep,col=r12))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Population size when the DMIs are resolved")+ylab("Prob(alive)")+labs(col = "Recombination rate")
ggplot(temp_data,aes(x=t_alive_both,y=p_alive/n_rep,col=log10(r12),pch=as.factor(arch)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the two DMIs")+ylab("Prob(alive)")+labs(col = "log(Recombination rate)") 
ggplot(temp_data,aes(x=p_alive/n_rep,y=(p_hybA1+p_hybB1)/n_rep,col=r12,pch=as.factor(arch)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+ylab("Prob(Hybrid speciation)")+xlab("Prob(alive)")+labs(col = "log(Recombination rate)")
ggplot(temp_data,aes(x=p_alive/n_rep,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=r12,pch=as.factor(arch)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+ylab("Prob(Hybrid speciation|alive)")+xlab("Prob(alive)")+labs(col = "log(Recombination rate)")
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=r12))+geom_point()+ xlab("Minimum population size")+ylab("Prob(alive)")+labs(col = "Recombination")+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom")+scale_colour_gradient(low = "orange", high = "purple")
Here we use a more precise estimation of the different outcomes
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==5000&whole_dataset_extra$arch!="sin_DMI",]Relationship between hybrid speciation, survival probability and recombination, with higher resolution
#pdf("Figure/cod_case_ep0d2_N5000_10k.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob. (Survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylim(0.4,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("Prob. (Hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.25) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob. (Hybrid sp.|Survival)")+xlab("Recombination")+geom_line()#+ guides(col=FALSE)
#dev.off()Relationship between minimum population size and hybrid speciation with higher resolution
#pdf("Figure/Nminvsrec_and_hyybvsrec_cod0d2_10k.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10(breaks=c(5*10^-4, 10^-2,0.5))+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))#+ guides(col=FALSE)
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("Prob.(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label)) + theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))#+ guides(col=FALSE)
#dev.off()Initial population size N=500
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==500,]#pdf("Figure/cod_case_ep0d2_N500.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylim(0,1)+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.15) + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()
#dev.off()ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
Initial population size of N=100
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$K==1000000&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$s1==-0.001&whole_dataset_extra$s2==-0.001&whole_dataset_extra$s3==-0.001&whole_dataset_extra$s4==-0.001,]#pdf("Figure/cod_case_ep0d2_N100.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob.(Survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob.(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+geom_line()## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()Initial population size of N=50000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==50000,]Relation between survival probability and hybrid speciation in large population
#pdf("Figure/cod_case_ep0d2_N50000.pdf",widt=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylim(0,1)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")
#dev.off()Juxtaposition of hybrid speciation and survival probability
For the Adjacent ABAB architecutre
temp_data=whole_dataset[whole_dataset$arch=="adj_ABAB"&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2)&whole_dataset$h12==1&whole_dataset$N %in% c(500,5000,50000),]#pdf("Figure/combined_surv_hyb_sp.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(N)))+geom_point()+geom_point(aes(x=r12,y=p_alive/n_rep,col=as.factor(N)),shape=0)+xlab("Recombination")+ylab("P(hybrid sp.) or P(alive)")+ scale_colour_manual( values=seq(8), name="Initial pop size") +theme(text = element_text(size=14))+ylim(0,1)+scale_x_log10()+theme(legend.position="bottom")
#ggplot(temp_data,aes(x=p_alive/n_rep,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch),labeller = as_labeller(arch_label)))+geom_point()+xlab("P(alive)")+ylab("P(hybrid sp.)")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylim(0,1)+ guides(col=FALSE)
#dev.off()###Investigating the effect of direct selection (or its lack) Same as above, but we are removing the effect of direct selection on all four loci
temp_data=whole_dataset_extra[whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$K==1000000&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$s1==0&whole_dataset_extra$s2==0&whole_dataset_extra$s3==0&whole_dataset_extra$s4==0,]#pdf("Figure/cod_case_ep0d2_N100_10k_nosel.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob.(Survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()+labs(tag = "B")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob.(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+geom_line()## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()temp_data=whole_dataset_extra[whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$K==1000000,]Comparaison between the case with and without direct selection
#pdf("Figure/comaprison_no_sel.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(s1)))+geom_point()+scale_x_log10()+ylab("Prob.(Survival)")+xlab("Recombination") +theme(text = element_text(size=14))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="right")+labs(col = "Selection")+labs(tags="A")#+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(s1)))+geom_point()+scale_x_log10()+ylab("Prob.(hybrid speciation)")+xlab("Recombination") +theme(text = element_text(size=14))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(s1)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation&survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(s1)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Selection")## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()Investigating the effect of K in small populations
temp_data=whole_dataset[whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==100&whole_dataset$arch!="sin_DMI"&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s3==-0.001&whole_dataset$s4==-0.001,]Same as above, but we are removing the varying the carrying capacity of the environment
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination") +theme(text = element_text(size=14))+facet_wrap(~arch)#+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination") +theme(text = element_text(size=14))+facet_wrap(~arch)#+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=14))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")+facet_wrap(~arch)+geom_smooth(se=F)+ guides(col=FALSE)## Warning in geom_point(se = F): Ignoring unknown parameters: `se`## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+facet_wrap(~arch)+geom_smooth(se=F)+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture")## Warning in geom_point(se = F): Ignoring unknown parameters: `se`## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))
# ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))
# ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_point()
# ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")
# ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_point(se=F)+scale_x_log10()+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation&survival)")+xlab("Recombination")+geom_smooth(se=F)Same, but with higher resolution
temp_data=whole_dataset_extra[whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$arch!="sin_DMI"&whole_dataset_extra$s1==-0.001&whole_dataset_extra$s2==-0.001&whole_dataset_extra$s3==-0.001&whole_dataset_extra$s4==-0.001,]#pdf("Figure/as_a_fucntion_of_K.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("Prob.(Survival)")+xlab("Recombination") +theme(text = element_text(size=14))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")#+ guides(col=FALSE)#+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture")
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(K)))+geom_point()+scale_x_log10()+ylab("Prob.(hybrid speciation)")+xlab("Recombination") +theme(text = element_text(size=14))+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation & survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(K)))+geom_point(se=F)+scale_x_log10() +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+facet_wrap(~arch,labeller =as_labeller( arch_label))+theme(legend.position="bottom")+labs(col = "Max pop. size")## Warning in geom_point(se = F): Ignoring unknown parameters: `se`
#dev.off()Effect of the mean offspring number
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==100000&whole_dataset$r12==0.1&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s3==-0.001&whole_dataset$s4==-0.001,]Effect of the mean offspring number on survival and hybrid specation in alrge populations
ggplot(temp_data,aes(x=mean_off,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Mean off")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+xlim(1,1.03)
ggplot(temp_data,aes(x=mean_off,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(survival)")+xlab("Mean off")+xlim(1,1.015)+geom_line()## Warning: Removed 89 rows containing missing values (`geom_point()`).## Warning: Removed 89 rows containing missing values (`geom_line()`).
ggplot(temp_data,aes(x=mean_off,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Mean off")+xlim(1,1.03)## Warning: Removed 5 rows containing missing values (`geom_point()`).
Effect of the mean offspring number on survival and hybrid specation in large populations
temp_data=whole_dataset[(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$N==1000&whole_dataset$r12==0.1&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s2==-0.001&whole_dataset$s4==-0.001,]ggplot(temp_data,aes(x=mean_off,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Mean off")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))
ggplot(temp_data,aes(x=mean_off,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(survival)")+xlab("Mean off")
ggplot(temp_data,aes(x=mean_off,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Mean off")## Warning: Removed 5 rows containing missing values (`geom_point()`).
Explore the effect of the population size N for each architecture
temp_data=whole_dataset[whole_dataset$mean_off==1.01&whole_dataset$K==1000000&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$r12==0.1&whole_dataset$s1==-0.001&whole_dataset$s2==-0.001&whole_dataset$s2==-0.001&whole_dataset$s4==-0.001,]Relationship between survival and hybrid speciation and the initial population size
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(survival)")+xlab("Initial hybrid population size")+scale_x_log10()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Initial hybrid population size")+scale_x_log10()
Zoom on the lower values for survival
#pdf("Figure/surv_vs_Nin_cod_ep0d2_r0d1.pdf",width=8,height=6)
ggplot(temp_data[temp_data$N<=1000,],aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob.(Survival)")+xlab("Initial hybrid population size")+scale_x_log10()+geom_line()
#dev.off()With higher resolution
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$s1==-0.001&whole_dataset_extra$s2==-0.001&whole_dataset_extra$s3==-0.001&whole_dataset_extra$s4==-0.001&whole_dataset_extra$K==1000000,]
temp_data$err_hyb=sqrt(((temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$p_alive*(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$p_alive)/temp_data$p_alive)#pdf("Figure/hybrid_spec_Nini.pdf",width=8, height=6)
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_line() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob. (Hybrid speciation|Survival)")+xlab("Initial size of the hybrid population")+scale_x_log10()+geom_point()+geom_errorbar(aes(x=N,ymin=(p_alive_hybA1+p_alive_hybB1)/p_alive-err_hyb,ymax=(p_alive_hybA1+p_alive_hybB1)/p_alive+err_hyb))
#dev.off()Decomposition of the outcome as a function of the initial population size
ggplot(temp_data,aes(x=N,y=(p_alive)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(alive)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+ylim(0.05,.2)+scale_x_log10(limits=c(10,750))## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 74 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 74 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation)")+xlab("Initial hybrid population size")+ scale_colour_manual(values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation)")+xlab("Initial hybrid population size")+ scale_colour_manual(values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_point()+ylim(0,0.05)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 35 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 35 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation & surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+ylim(0,0.05)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 10 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 10 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A1B2 & surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybB1)/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1 & surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_alive_both,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1 & surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_alive_both,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Time resolution of DMIs")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+scale_y_log10()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_diff,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_diff_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs|alive")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=t_diff_alive_hybA1,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs|alive for hyb A1")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 1 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=t_diff_alive_hybB1,col=as.factor(arch)))+geom_smooth(se=F)+ylab("Difference in resolution of DMIs|alive for hyb B1")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(Alive)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid sp.|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_line()+ylab("P(hybrid sp.|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A1B2|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+ylim(0,0.05)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 1 rows containing missing values (`geom_smooth()`).
ggplot(temp_data,aes(x=N,y=(p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+ylim(0,0.05)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning: Removed 20 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 20 rows containing missing values (`geom_point()`).
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+ylab("P(Alive)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+geom_line()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+ylab("P(hybrid speciation|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()+geom_smooth(se=F)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A1B2|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(temp_data,aes(x=N,y=(p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_smooth(se=F)+ylab("P(hybrid speciation A2B1|surviving)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_point()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Stacks plot representing the different outcomes for the different
architectures
temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="adj_ABBA"&!is.na(whole_dataset_extra$arch)&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_adjABBA.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Adjacent ABBA")
#dev.off()temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="adj_ABAB"&!is.na(whole_dataset_extra$arch)&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_adjABAB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Adjacent ABAB")
#dev.off()temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="cro_AABB"&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_croAABB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Crossed AABB")
#dev.off()temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="cro_ABBA"&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_croABBA.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Crossed ABBA")
#dev.off()temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="nes_AABB"&whole_dataset_extra$K==10^6&whole_dataset_extra$n_rep==10000&whole_dataset_extra$s1==-0.001,]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_nesAABB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Nested AABB")
#dev.off()temp_data=whole_dataset_extra[(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01&whole_dataset_extra$arch=="nes_ABAB",]
temp_data=temp_data[order(temp_data$N),]
#pdf("Figure/stacked_plot_nesABAB.pdf",width=8,height=6)
stack_plot(temp_data,ymax = 1,title = "Nested ABAB")
#dev.off()Decomposition with higher resolution
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==1&whole_dataset_extra$r12==0.1&whole_dataset_extra$mean_off==1.01,]
ggplot(temp_data,aes(y=(p_alive_hybA1+p_alive_hybB1)/n_rep,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_y_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(Hybrid speciation&survival)")+xlab("Population size at the time of resolution of the genetic conflicts")
#ggplot(temp_data,aes(y=NL_alive_A1,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_y_log10()
#ggplot(temp_data,aes(y=NL_alive_A2,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_y_log10()
#ggplot(temp_data,aes(y=NL_alive_B1,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_y_log10()
#ggplot(temp_data,aes(y=NL_alive_B2,x=N_alive_both,col=arch))+geom_point()+scale_x_log10()+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_y_log10()#temp_soft=read.table("~/IGC/dmi_and_pattern/temp_soft_sel.csv",sep=",",header=TRUE)
#ggplot(temp_soft,aes(x=N,y=(p_hybA1 +p_hybB1)/1000,col=arch ))+geom_point()+theme(text = element_text(size=14))+scale_x_log10()+ylim(0,0.2)+ylab("Hybdrid spe.")+geom_line()Relation between survival and initial population size for unlinked loci
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$r12==0.5&whole_dataset$mean_off==1.01,]ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(survival)")+xlab("Initial hybrid population size")+scale_x_log10()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Initial hybrid population size")+scale_x_log10()
Relation between survival and initial population size for closely linked loci
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==1&whole_dataset$r12==0.01&whole_dataset$mean_off==1.01,]ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("P(hybrid speciation)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+scale_x_log10()
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(survival)")+xlab("Initial hybrid population size")+scale_x_log10()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Initial hybrid population size")+scale_x_log10()
Lethal epistasis: epistasis is fixed to e=-0.99
Initial population size of N=5000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==1&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]Relation between survival, hybrid speciation and recombination for (almost) lethal epistasis
#pdf("Figure/cod_case_ep0d99_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=(p_alive)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("Prob.(survival)")+xlab("Recombination rate")+scale_colour_manual( values=seq(7), labels =list_Arch_legend ,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line() #+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("Prob.(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.25) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+ guides(col=FALSE)
#dev.off()Relation between survival and minimum population size for (almost) lethal epistasis
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,0.2)+ylab("P(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))
####N=50000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==1&whole_dataset$N==50000,]Relation between survival, hybrid speciation and recombination for (almost) lethal epistasis in large populations
ggplot(temp_data,aes(x=r12,y=(p_alive)/n_rep,col=as.factor(arch)))+geom_point()+geom_smooth()+scale_x_log10()+ylim(0,1)+ylab("Survival probability")+xlab("Recombination rate")+scale_colour_manual( values=seq(7), labels = list_Arch,name="Linkage\n architecture") ## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
N=100
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.99|whole_dataset_extra$e13==-0.99|whole_dataset_extra$e14==-0.99)&whole_dataset_extra$h12==1&whole_dataset_extra$N==100&whole_dataset_extra$arch!="sin_DMI",]Relation between survival, hybrid speciation and recombination for (almost) lethal epistasis in small populations with higher resolution
#pdf("Figure/cod_case_ep0d99_N100.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))
#dev.off()#Analysis of the linked loci: Recessive scenario
All architectures
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&whole_dataset$h12==0&whole_dataset$arch!="sin_DMI",]Relationship between survival probability and minimum population size
#ggplot(temp_data,aes(x=t_alive_both,y=p_alive/n_rep,col=log10(r12),pch=as.factor(arch),labeller = as_labeller(arch_label)))+geom_point()+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Time of resolution of the two DMIs")+ylab("Prob(alive)")+labs(col = "log(Recombination rate)")+scale_x_log10()
ggplot(temp_data,aes(x=N_min_alive,y=p_alive/n_rep,col=r12))+geom_point(size=0.5)+ scale_colour_gradient(low = "orange", high = "purple")+xlab("Minimum population size")+ylab("P(survival)")+labs(col = "Recombination rate")+scale_x_log10(breaks=c(10,1000,100000))+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(legend.position="bottom",legend.key.size = unit(0.23, "in"))+geom_segment(aes(x = 10, y = 0, xend = 600, yend = 1),size=1,color="black")
Epistasis is fixed to e=-0.2
temp_data=whole_dataset[default_regime(whole_dataset)&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==0&whole_dataset$r12==0.1&whole_dataset$mean_off==1.01,]
temp_data$err_hyb=sqrt(((temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$p_alive*(temp_data$p_hybA1+temp_data$p_hybB1)/temp_data$p_alive)/temp_data$p_alive)Relationship between survival, hybrid speciation and initial population size for all 6 architectures
#pdf("Figure/rec_case_all_arch_vs_Nini.pdf",width=8,height=6)
ggplot(temp_data,aes(x=N,y=p_alive/n_rep,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob.(Survival)")+xlab("Initial hybrid population size")+scale_x_log10()+geom_line()
ggplot(temp_data,aes(x=N,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+ylab("Prob.(hybrid speciation)")+xlab("Initial hybrid population size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+scale_x_log10()+geom_line()
ggplot(temp_data,aes(x=N,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Initial hybrid population size")+scale_x_log10()+geom_line()+geom_errorbar(aes(x=N,ymin=(p_alive_hybA1+p_alive_hybB1)/p_alive-err_hyb,ymax=(p_alive_hybA1+p_alive_hybB1)/p_alive+err_hyb))
#dev.off()Relationship between survival, hybrid speciation and recombination for all 6 architectures
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.2|whole_dataset$e13==-0.2|whole_dataset$e14==-0.2)&whole_dataset$h12==0&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]#pdf("Figure/rec_case_ep0d2_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(alive)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylim(0.4,1)+geom_line()# guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()# guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()# guides(col=FALSE)
#dev.off()Same as above with higher resolution
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$h12==0&whole_dataset_extra$N==5000&whole_dataset_extra$arch!="sin_DMI",]#pdf("Figure/rec_case_ep0d2_N5000_10k.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob.(Survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylim(0.85,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("Prob.(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()# guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+ geom_line()#guides(col=FALSE)
#dev.off()Relationship between recombination, minimum population size and hybrid speciation, with higher resolution
#pdf("Figure/Nminvsrec_and_hyybvsrec_rec0d2_10k.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=N_min_alive,col=as.factor(arch)))+geom_point()+scale_x_log10(breaks=c(5*10^-4, 10^-2,0.5))+ylab("Min. pop. size")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label))+ theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))
ggplot(temp_data,aes(x=N_min_alive,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10(breaks=c(500,1000, 2000,4000))+ylim(0,0.65)+ylab("Prob.(hybrid speciation)")+xlab("Min. pop. size")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ guides(col=FALSE)+facet_wrap(~as.factor(arch),labeller = as_labeller(arch_label)) + theme(axis.text.x = element_text(angle=45),plot.margin=margin(5.5,15,0,0,"pt"))
#dev.off()Recessive and codominant together
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$N==5000&whole_dataset_extra$arch!="sin_DMI",]Relationship between hybrid speciation, survival probability ansd recombination, with higher resolution
#pdf("Figure/case_ep0d2_N5000_10k.pdf",width=8,height=6)
temp_data$dom=""
temp_data$dom[temp_data$h12==1]="Codominant"
temp_data$dom[temp_data$h12==0]="Recessive"
temp_data$extra=plyr::mapvalues(temp_data$dom, c("Codominant", "Recessive"), c("A","B"))
temp_data$extra2=plyr::mapvalues(temp_data$dom, c("Codominant", "Recessive"), c("A",""))
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+geom_text(x=-4,y=0.75,aes(label=extra),colour="black")+geom_text(x=-4,y=0.99,aes(label=extra),colour="black")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()+facet_grid(rows=vars(dom),scales="free")#+theme(legend.position="bottom")#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob. (Hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()+facet_grid(rows=vars(dom),scales="free")+geom_text(x=-4,y=0.6,aes(label=extra),colour="black")+geom_text(x=-4,y=0.15,aes(label=extra2),colour="black")#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10() + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob. (Hybrid sp.|Survival)")+xlab("Recombination")+geom_line()+facet_grid(rows=vars(dom),scales="free")+geom_text(x=-4,y=0.6,aes(label=extra),colour="black")+geom_text(x=-4,y=0.22,aes(label=extra2),colour="black")#+ guides(col=FALSE)
#dev.off()#pdf("Figure/case_ep0d2_N5000_10k_v2.pdf",width=8,height=6)
temp_data=whole_dataset_extra[default_regime(whole_dataset_extra)&whole_dataset_extra$mean_off==1.01&(whole_dataset_extra$e12==-0.2|whole_dataset_extra$e13==-0.2|whole_dataset_extra$e14==-0.2)&whole_dataset_extra$N==5000&whole_dataset_extra$arch!="sin_DMI",]
temp_data2=temp_data
temp_data2$hybrid_spe=(temp_data2$p_hybA1+temp_data2$p_hybB1)/temp_data$n_rep
temp_data2$conditioned=0
temp_data3=temp_data
temp_data3$hybrid_spe=(temp_data3$p_alive_hybA1+temp_data3$p_alive_hybB1)/temp_data3$p_alive
temp_data3$conditioned=1
temp_data=rbind(temp_data2,temp_data3)
dom_lab=c(`0`="Recessive",`1`="Codominant")
names(dom_lab)=c("0","1")
conditioned_lab=c(`0`="Unconditional probability",`1`="Conditional probability")
names(conditioned_lab)=c("0","1")
temp_data$Cond=sapply(c(1:length(temp_data[,1])), function(x) paste(temp_data$h12[x], temp_data$conditioned[x], sep="_"))
temp_data$extra=plyr::mapvalues(temp_data$Cond, c("0_0", "0_1" , "1_0","1_1"), c("A","B","C","D"))
temp_data$extra2=plyr::mapvalues(temp_data$Cond, c("0_0", "0_1" , "1_0","1_1"), c("","","C","D"))
ggplot(temp_data,aes(x=r12,y=hybrid_spe,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob. (Hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+geom_line()+facet_grid(h12~conditioned,labeller =labeller(.rows=dom_lab,.cols=conditioned_lab),scales="free")+theme(legend.position="bottom")+geom_text(x=-4.3,y=0.675,aes(label=extra),colour="black")+geom_text(x=-4.3,y=0.225,aes(label=extra2),colour="black")+theme(text = element_text(size=14))+theme(panel.spacing.x = unit(1, "lines")) +
coord_cartesian(clip = "off")
#dev.off()Lethal epistasis: epistasis is fixed to e=-0.99
###N=5000
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==0&whole_dataset$N==5000&whole_dataset$arch!="sin_DMI",]Relationship between survival, hybrid speciation and recombination for all 6 architectures for lethal epistasis
#pdf("Figure/rec_case_ep0d99_N5000.pdf",width=8,height=6)
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("Prob.(Survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylim(0.85,1)+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("Prob.(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ geom_line()#guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\narchitecture") +theme(text = element_text(size=14))+ylab("Prob.(hybrid speciation|survival)")+xlab("Recombination")+geom_line()# guides(col=FALSE)
#dev.off()###N=500
temp_data=whole_dataset[default_regime(whole_dataset)&whole_dataset$mean_off==1.01&(whole_dataset$e12==-0.99|whole_dataset$e13==-0.99|whole_dataset$e14==-0.99)&whole_dataset$h12==0&whole_dataset$N==500&whole_dataset$arch!="sin_DMI",]Relationship between survival, hybrid speciation and recombination for all 6 architectures for lethal epistasis and small populatiobs
ggplot(temp_data,aes(x=r12,y=p_alive/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylab("P(survival)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+geom_line()#+ guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_hybA1+p_hybB1)/n_rep,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1)+ylab("P(hybrid speciation)")+xlab("Recombination")+ scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ geom_line()#guides(col=FALSE)
ggplot(temp_data,aes(x=r12,y=(p_alive_hybA1+p_alive_hybB1)/p_alive,col=as.factor(arch)))+geom_point()+scale_x_log10()+ylim(0,1) + scale_colour_manual( values=seq(7), labels = list_Arch_legend,name="Linkage\n architecture") +theme(text = element_text(size=14))+ylab("P(hybrid speciation|survival)")+xlab("Recombination")+geom_line()# guides(col=FALSE)
#Evolutionary rescue
Here we simulated with the same program an evolutionary rescue scenario, and evaluated the link between minimum population size and survival probability.
Reading the simulations output and writing the dataframe
if (FALSE){
fake_dataset=as.data.frame(rbind(rep(NA,length(names_col)),rep(NA,length(names_col))))
colnames(fake_dataset)=names_col
} else{
fake_dataset=read.csv("evol_rescue.csv")
fake_dataset$arch=as.factor(list_Arch[fake_dataset$arch])
}known_folder_f=read.table("know_folder_f.txt")
list_folder_f=list.dirs(path = ".", full.names = TRUE, recursive = TRUE)
list_folder_f=gsub("./","",list_folder_f)
list_folder_f=list_folder_f[grep("^fake", list_folder_f)]update=FALSE
#known_folder_f=c()
for (i in list_folder_f){
if (!(i %in% unlist(known_folder_f))){
update=TRUE
#check_non_finish(i)
fake_dataset=add_dataframe(i,fake_dataset)
known_folder_f=c(known_folder_f,i)
}
}
#whole_dataset=whole_dataset[!is.na(whole_dataset$n_rep),]
if (update){
write.table(list_folder_f,"know_folder_f.txt")
write.csv(fake_dataset,"evol_rescue.csv",row.names = FALSE)
} fake_dataset=read.csv("evol_rescue.csv")
fake_dataset$arch=as.factor(list_Arch[fake_dataset$arch])
fake_dataset=fake_dataset[-c(1,2),]
fake_dataset$N_min_alive[is.na(fake_dataset$N_min_alive)]=0##Analysis
Link between survival (or evolutionary rescue) and minimum population size
#pdf("Figure/prob_rescue.pdf",width=8,height=6)
ggplot(fake_dataset,aes(x=N_min_alive,y=p_alive/n_rep,col=as.factor(f1000/N)))+geom_point()+scale_x_log10(limits=c(1,10^5))+ geom_segment(aes(x = 10, y = 0, xend = 600, yend = 1),size=.5,color="black")+ylab("Prob.(Survival)")+xlab("Minimum population size")+labs(col="Freq. of\n allele A")+facet_wrap(~as.factor(signif(fake_dataset$s1,3)),nrow=3)## Warning: Transformation introduced infinite values in continuous x-axis
#dev.off()Link between strength of selection and minimum population size
ggplot(fake_dataset,aes(x=s1,y=N_min_alive,col=as.factor(f1000/N)))+geom_point()+ylab("Min pop size")+xlab("Selection on the ancestral allele")+labs(col="Freq. of allele A")+geom_line()
Deterministic trajectories
if(TRUE){
case_adjabab=read.table("deterministic_output/freq_sym_adjabab.txt")
case_adjabba=read.table("deterministic_output/freq_sym_adjabba.txt")
case_croaabb=read.table("deterministic_output/freq_sym_croaabb.txt")
case_croabba=read.table("deterministic_output/freq_sym_croabba.txt")
case_nesabab=read.table("deterministic_output/freq_sym_nesabab.txt")
case_nesaabb=read.table("deterministic_output/freq_sym_nesaabb.txt")
} else {
case_adjabab=read.table("deterministic_output/freq_asymA_adjabab.txt")
case_adjabba=read.table("deterministic_output/freq_asymA_adjabba.txt")
case_croaabb=read.table("deterministic_output/freq_asymA_croaabb.txt")
case_croabba=read.table("deterministic_output/freq_asymA_croabba.txt")
case_nesabab=read.table("deterministic_output/freq_asymA_nesabab.txt")
case_nesaabb=read.table("deterministic_output/freq_asymA_nesaabb.txt")
}case_adjabab=cbind(case_adjabab,rep("Adj. ABAB",1001))
case_adjabba=cbind(case_adjabba,rep("Adj. ABBA",1001))
case_croaabb=cbind(case_croaabb,rep("Cro. AABB",1001))
case_croabba=cbind(case_croabba,rep("Cro. ABAB",1001))
case_nesabab=cbind(case_nesabab,rep("Nes. ABAB",1001))
case_nesaabb=cbind(case_nesaabb,rep("Nes. AABB",1001))colnames(case_adjabab)=c("pop_size","freq_A1","freq_B1","freq_A2","freq_B2","freq_anc","freq_ancB2","freq_ancA2","freq_D2","freq_ancB1","freq_parB","freq_hybB1","freq_D2B","freq_ancA1","freq_hybA1","freq_parA","freq_D2A","freq_D1","freq_D1B","freq_D1A","freq_DD")
colnames(case_croaabb)=c("pop_size","freq_A1","freq_A2","freq_B1","freq_B2","freq_anc","freq_ancB2","freq_ancB1","freq_parB","freq_ancA2","freq_D2","freq_hybB1","freq_D2B","freq_ancA1","freq_hybA1","freq_D1","freq_D1B","freq_parA","freq_D2A","freq_D1A","freq_DD")
colnames(case_nesabab)=c("pop_size","freq_A1","freq_B2","freq_A2","freq_B1","freq_anc","freq_ancB1","freq_ancA2","freq_hybB1","freq_ancB2","freq_parB","freq_D2","freq_D2B","freq_ancA1","freq_D1","freq_parA","freq_D1A","freq_hybA1","freq_D1B","freq_D2A","freq_DD")
colnames(case_adjabba)=c("pop_size","freq_A1","freq_B1","freq_B2","freq_A2","freq_anc","freq_ancA2","freq_ancB2","freq_D2","freq_ancB1","freq_hybB1","freq_parB","freq_D2B","freq_ancA1","freq_parA","freq_hybA1","freq_D2A","freq_D1","freq_D1A","freq_D1B","freq_DD")plot(-100,xlim=c(0,500),ylim=c(0,500),ylab="Expected population size",xlab="Generations")
points(seq(0,1000),case_adjabab$V1,col="black",type="l")
points(seq(0,1000),case_adjabba$V1,col="red",type="l")
points(seq(0,1000),case_croaabb$V1,col="green",type="l")
points(seq(0,1000),case_croabba$V1,col="blue",type="l")
points(seq(0,1000),case_nesabab$V1,col="purple",type="l")
points(seq(0,1000),case_nesaabb$V1,col="cyan",type="l")
plot(-1,xlim=c(0,50),ylim=c(0,0.7),ylab="Freq. ancestral hap.",xlab="Generations")
points(seq(0,1000),case_adjabab$V6,col="black",pch=16)
points(seq(0,1000),case_adjabba$V6,col="red",pch=16)
points(seq(0,1000),case_croaabb$V6,col="green",pch=16)
points(seq(0,1000),case_croabba$V6,col="blue",pch=16)
points(seq(0,1000),case_nesabab$V6,col="purple",pch=16)
points(seq(0,1000),case_nesaabb$V6,col="cyan",pch=16)
plot(NA,xlim=c(0,.7),ylim=c(50,5000),xlab="Freq. ancestral hap.",ylab="Expected population size",log="y")
points(case_adjabab$V6[1:51],case_adjabab$V1[1:51],col="black",pch=16)
points(case_adjabba$V6[1:51],case_adjabba$V1[1:51],col="red",pch=16)
points(case_croaabb$V6[1:51],case_croaabb$V1[1:51],col="green",pch=16)
points(case_croabba$V6[1:51],case_croabba$V1[1:51],col="blue",pch=16)
points(case_nesabab$V6[1:51],case_nesabab$V1[1:51],col="purple",pch=16)
points(case_nesaabb$V6[1:51],case_nesaabb$V1[1:51],col="cyan",pch=16)
Below we are representing the frequencies of haplotypes or group of haplotypes using the following color scheme, black for ancestral (a1b1a2b2), blue (a1B1a2B2) and red (A1b1A2b2) for the parental haplotypes, pink (A1b1a2b2 or a1b1A2b2) and cyan (a1B1a2b2 or a1b1a2B2), green (a1B1A2b2) and yellow (A1b1a2B2) for the hybrid speciation haplotypes and purple for any haplotype harboring at least one DMI.
temp_case=case_adjabab
time_step=50
plot(-100,xlim=c(0,time_step),ylim=c(0,1),ylab="Freq.(hap).",xlab="Generations",cex.lab=1.5,cex.axis=1.5, main ="Adjacent ABAB")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(rep(0,(time_step+1)),temp_case$freq_anc[seq((time_step+1),1,-1)]),col="black")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(temp_case$freq_anc[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA)[seq((time_step+1),1,-1)]),col="red")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[seq((time_step+1),1,-1)]),col="blue")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[seq((time_step+1),1,-1)]),col="pink")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[seq((time_step+1),1,-1)]),col="cyan")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[seq((time_step+1),1,-1)]),col="yellow")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[seq((time_step+1),1,-1)]),col="green")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[1:(time_step+1)],(rep(1,time_step+1))[seq((time_step+1),1,-1)]),col="purple")
temp_case=case_croaabb
time_step=50
plot(-100,xlim=c(0,time_step),ylim=c(0,1),ylab="Freq.(hap).",xlab="Generations",cex.lab=1.5,cex.axis=1.5,main="Crossed AABB")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(rep(0,(time_step+1)),temp_case$freq_anc[seq((time_step+1),1,-1)]),col="black")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(temp_case$freq_anc[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA)[seq((time_step+1),1,-1)]),col="red")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[seq((time_step+1),1,-1)]),col="blue")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[seq((time_step+1),1,-1)]),col="pink")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[seq((time_step+1),1,-1)]),col="cyan")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[seq((time_step+1),1,-1)]),col="yellow")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[seq((time_step+1),1,-1)]),col="green")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[1:(time_step+1)],(rep(1,time_step+1))[seq((time_step+1),1,-1)]),col="purple")
temp_case=case_nesabab
time_step=50
plot(-100,xlim=c(0,time_step),ylim=c(0,1),ylab="Freq.(hap).",xlab="Generations",cex.lab=1.5,cex.axis=1.5,main="Nested ABAB")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(rep(0,(time_step+1)),temp_case$freq_anc[seq((time_step+1),1,-1)]),col="black")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(temp_case$freq_anc[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA)[seq((time_step+1),1,-1)]),col="red")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[seq((time_step+1),1,-1)]),col="blue")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[seq((time_step+1),1,-1)]),col="pink")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[seq((time_step+1),1,-1)]),col="cyan")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[seq((time_step+1),1,-1)]),col="yellow")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[seq((time_step+1),1,-1)]),col="green")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[1:(time_step+1)],(rep(1,time_step+1))[seq((time_step+1),1,-1)]),col="purple")
temp_case=case_adjabba
time_step=50
plot(-100,xlim=c(0,time_step),ylim=c(0,1),ylab="Freq.(hap).",xlab="Generations",cex.lab=1.5,cex.axis=1.5,main="Adjacent ABBA")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(rep(0,(time_step+1)),temp_case$freq_anc[seq((time_step+1),1,-1)]),col="black")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c(temp_case$freq_anc[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA)[seq((time_step+1),1,-1)]),col="red")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[seq((time_step+1),1,-1)]),col="blue")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[seq((time_step+1),1,-1)]),col="pink")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[seq((time_step+1),1,-1)]),col="cyan")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[seq((time_step+1),1,-1)]),col="yellow")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1)[1:(time_step+1)],(temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[seq((time_step+1),1,-1)]),col="green")
polygon(c(seq(0,time_step),seq(time_step,0,-1)),c((temp_case$freq_anc+temp_case$freq_parA+temp_case$freq_parB+temp_case$freq_ancA1+temp_case$freq_ancA2+temp_case$freq_ancB1+temp_case$freq_ancB2+temp_case$freq_hybA1+temp_case$freq_hybB1)[1:(time_step+1)],(rep(1,time_step+1))[seq((time_step+1),1,-1)]),col="purple")